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ABSTRACT: We investigate the dynamics of the collapse of a single copolymer chain, when the 
solvent quality is suddenly quenched from good to poor. We employ Brownian dynamics simu- 
lations of a bead-spring chain model and incorporate fluctuating hydrodynamic interactions via 
the Rotne-Prager-Yamakawa tensor. Various copolymer architectures are studied within the 
framework of a two-letter HP model, where monomers of type H (hydrophobic) attract each 
other, while all interactions involving P (polar or hydrophilic) monomers are purely repulsive. 
The hydrodynamic interactions are found to assist the collapse. Furthermore, the chain se- 
quence has a strong influence on the kinetics and on the compactness and energy of the final 
state. The dynamics is typically characterised by initial rapid cluster formation, followed by 
coalescence and final rearrangement to form the compact globule. The coalescence stage takes 
most of the collapse time, and its duration is particularly sensitive to the details of the archi- 
tecture. Long blocks of type P are identified as the main bottlenecks to find the globular state 
rapidly. 

1. Introduction 

The dynamics of the conformational change of a single polymer chain from a swollen coil to a collapsed 
globule due to a sudden change in the solvent quality has received much attention over the past 
decadesJMim 0ne of the primary reasons for the continued interest in polymer collapse dynamics is 
because of its qualitative similarity with the folding transitions seen in proteinsP It is believed that 
understanding the collapse transition will enable us to obtain better insight into the conformational 
transitions in many complex biological systems, and to understand how other bio-molecules such as 
DNA react to a change in their environment. 
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Since the folding or collapse process of a chain mainly occurs in an environment surrounded by 
solvent molecules, it is now widely recognized that it is important to take into account the effects of the 
solvent-mediated long-range dynamic correlations between different segments of the chain, known as 
hydrodynamic interactions (HI). Although HI does not affect static properties, it seriously alters the 
dynamic properties of semi-dilute and dilute solutions. For instance, a recent study^ of the diffusion 
and folding of several model proteins has shown that molecular simulations which neglect HI are 
incapable of reproducing the expected experimental translational and rotational diffusion coefficients 
of folded proteins, while simulations that include HI predict the expected experimental values very 
well. The study of the folding process of small proteins has also revealed that inclusion of HI hastens 
the folding process by at least a factor of twoP^ The analogous and simpler problem of the kinetics of 
homopolymer collapse both in the absence and presence of HI has been studied via a variety of different 
theoretical and numerical approachesl^al 9 ' 10 ' 11 ' 12 ' 4 ' 13 ' 14 ' 15 ' 10 ' 17 ' 18 ' 19 ' 20 ' ^ 1 ' ^ ' ^ ' 24 ' and it is probably fair 
to say that a reasonable degree of understanding has been obtained. 

Although the homopolymer model has been widely used as a prototype to understand the protein 
folding problem, it cannot capture all its complexities. For instance, an important missing detail in a 
homopolymer model is the formation of hydrophobic residues (H) or amino acids in the core and polar 
residues (P) on the outer surface of a collapsed globule. Moreover, it is believed that the kinetic ability 
of natural proteins to fold results from the evolutionary selection of the sequence of the molecule's 
amino acids, which should also result in a unique native structure, i. e. a collapsed state of vanishing 
or at least low degeneracy in terms of its (free) energy. These are aspects that cannot be captured 
by a homopolymer model; rather at least two different types of interactions or "letters" are needed to 
provide the possibility to encode a native structure. 

For convenience as well as efficiency in numerical calculations, we choose to use a simple two-letter 
code HP model to study the kinetics of collapse of heteropolymers. Such models have been referred to 
as prototeinsf^ In this model, a chain is represented by a binary sequence of H (hydrophobic) and P 
(polar) monomers. Although the HP model has been widely used to study the protein folding problem 
as well as heteropolymers in general, almost all the results are obtained from predictions on discrete 
lattice models, ^ 5 * 26 * 27 ! while direct simulation studies in the continuum are scarce.^ From such lattice 
models one can conclude that only a small fraction of sequences folds uniquely. Furthermore, not all 
geometric native structures can be encoded in terms of a suitably adjusted sequence! 26 ! 27 ! 

For chains of reasonable length, the number of possible sequences and conformations is too large for 
exact enumeration and hence statistical sampling is required. In the present study, we restrict attention 
to sequences that contain 50% H monomers and 50% P monomers; there are both biologicaP^ as well as 
theoretica l 30 ! 31 ! arguments that this case should be particularly relevant for protein folding. Previous 
theoretical and simulational studies of heteropolymer collapsj 32 ! 33 ! 34 ! 35 ! 36 ! 37 ! 38 ! 2 ^ have shown that it 
typically comprises two stages, namely the rapid formation of clusters or micelles along the chain, 
followed by coalescence, whose speed is rather sensitive to the details of the sequence. Furthermore, 
the final geometry depends on these details as well and is not always a compact spherical globule. 

With the hope of gaining insight into the protein folding problem, there has been a trend in recent 
years to design or engineer (in theoretical or computer models) sequences that mimic biological evo- 
lution such that they rapidly fold to a desired target native structure! 39 ! 40 ! 41 ! In particular, Khokhlov 
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and Khalatui * 36 * 37 ! have introduced such a method for a copolymer chain composed of 50% H and 50% 
P monomers. It involves first collapsing a homopolymer chain into its equilibrium compact state and 
then marking half of the monomers that are closest to the centre of mass as H type while the remaining 
monomers are marked as P type. Copolymer chains with sequences generated via this method (also 
known as "protein-like" copolymers) were found to fold much more rapidly and the native confor- 
mations were more stable compared to random copolymers and random block copolymers with the 
same average block length of H or P polymers. This result was however obtained for a lattice model 
without inclusion of HI. 

In the present study, we investigate the same question for an off-lattice bead-spring model under 
both absence and presence of HI. Furthermore, we investigate the influence of the length of H and P 
blocks, i. e. we compare the "protein-like" copolymers (PLC) not only to random block copolymers 
with the same average block length (RBC), but also to multi-block copolymers (MBC) with various 
different block lengths. Finally, wc attempt to identify the key features of an HP sequence that 
govern the collapse kinetics. In order to capture realistic dynamics, we use Brownian Dynamics (BD) 
where HI is taken into account via the Rotne-Prager-Yamakawa (RPY) tensor. To the best of our 
knowledge, the combined effects of HI and chain sequence on collapse kinetics within the framework 
of BD simulations has not been reported so far in the literature. 

In Sec. [21 we provide the basic equations for the bead-spring chain model and the BD simulation 
method. Section [3] discusses our key findings, which are summarised in Sec. [4] 

2. Simulation method 

2.1. Model and equation of motion 

The macromolecule is represented by a bead-spring chain model of N beads, which are connected by 
N — 1 springs. Its configuration is specified by the set of position vectors r M (/x = 1,2,.. .N). The 
dynamics is governed by the Ito stochastic differential equation^ 

r;(t* + At*) = r* (f) + jD^ • (f s / + Ff*) At* + -±=B^ • AW„ v — 1,2, ... N. (1) 

Here we use dimensionless units with Ih = yk^T/H and A# = (^/AH as elementary length scale 
and time scale, respectively w here fcs is the Boltzmann constant, T the temperature, H the 
spring constant and C (C = 67rr? s a, where r] s is the solvent viscosity) the Stokes friction coefficient of 
a spherical bead of radius a. Dimensionless quantities are introduced via r* = t^/Ih and t* = t/Xn- 
W„ is a Wiener process, D M „ is the diffusion tensor representing the effect of the motion of a bead fj, 
on another bead v and is defined as D M „ = S^ u 5 + fi^, where 5^ is the Kronecker delta, S is the unit 
tensor and fi M „ is the hydrodynamic interaction tensor. F^,* is the dimensionless spring force, F" 1 '* 
is the dimensionless excluded-volume force and the components of B M „ are related to the HI tensor 
such that D M „ = -Bj^P^l Throughout, the summation convention is implied for repeated indices. 
We use the regularised Rotne-Prager-Yamakawa (RPY) tensor to represent Hlj 45 * 46 ! its form is 
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where ft* is the dimensionless bead radius in the bead-spring model, defined as ft* = a/(///y / 7r). 
Typical values of ft* lie between and 0.5. 

We have employed a two-letter code HP model for our copolymers, each chain being composed of 
hydrophobic H and polar P beads. The polymer is assumed to be in an aqueous solvent such that HP 
and PP interactions are purely repulsive (i. e. in good solvent), while HH interactions are attractive 
(i. e. in poor solvent). Similar to previous workj^we use the dimensionless potentials Vev for purely 
repulsive and V a ttr for attractive interactions; they are defined as 
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Here a and £lj are the Lennard- Jones parameters: e* = (-u/k^T is the dimensionless strength of 
the excluded volume interaction or quench depth, which was used as a variable parameter; in most 
studies, we chose the value 2.5 (which is always assumed unless otherwise stated). The parameter 
ct/Ih denotes the corresponding dimensionless length scale, which we set at \fl . The dimensionless 
distance from bead [i to bead v is denoted by r*„. The cutoff radius of Kittr is R* — 2.5 (ct/Ih) 
and the function c(R* c ) is chosen such that the value of the potential is zero at the cutoff, i.e., 
c(R* c ) = [(a/l H )/Rtr-[(a/l H )/R* c } e . 

The adjacent beads in the chain also interact via a finitely extensible non- linear elastic (FENE) 
spring potential 

^FENE(r*) = -ilnfl-^V (8) 



where we choose Qq = 2 (ct/Ih) = 2^/7 for the dimensionless maximum stretchable length of a single 
spring. 

The strength of HI is governed by the dimensionless Stokes radius h*. Motivated by previous 
work on FENE chains J^ 1 ^ 1 we choose h* = 0.5^/28/33 w 0.46; this choice defines the HI strength as 
h* = 0.5 in units of the FENE equilibrium bond lengthPSM Our main guidance in choosing h* = 0.5 
came from previous results on homopolymers™ where this value led to a rapid collapse. Similarly, 
we found in Ref. [M] that for e^/k-oT = 2.5 the homopolymer chain smoothly folds into its final 
equilibrium compact stage without being trapped in one of its intermediate mctastable states. 

For all simulations a time step size of At* = 0.001 was used, yielding a solution with very small 
discretisation errors. For further technical details on the BD algorithm, see Ref. [43j 
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The simulation procedure was carried out as a "quenching" experiment where the chain was started 
in a conformation corresponding to good solvent conditions. The starting conformation was obtained 
by equilibrating a homopolymer where all interactions are purely repulsive. This equilibration was 
done via a BD run without HI for a period of T cq = 15 t^ r , where r^p = 0.5 sin _2 (7r/2iV) is the longest 
dimcnsionless Rouse relaxation time, starting from a random walk configuration. This procedure 
was repeated roughly 500 times in order to generate a sample of statistically independent starting 
configurations which were stored for later use in the quenching runs. At time t = 0, a sequence of H 
and P blocks was selected (see below), and attractive interactions between HH pairs were turned on. 
From then on, the collapse process was monitored. 

2.2. Chain sequence construction 

We have carried out simulations for two different chain lengths of N = 64 and 128 beads at a fixed 
H:P ratio of 1 : 1 (i. c. Ah = Np = N/2). We have studied three different families of copolymer 
chains with sequence types that were either inherited from the parent globule, regular or probabilistic. 

The first sequence type is the protein-like copolymer or PLC, where the sequence has been gener- 
ated by the following process. A homopolymer is picked from the set of 500 stored initial conforma- 
tions, and then collapsed to a compact globule by making all the beads attractive. This is followed 
by marking the interior as hydrophobic H monomers and the exterior as hydrophilic or polar P ones. 
Typically, this procedure leads to a sequence of H and P blocks which look at least partially disor- 
dered when viewed along the chain backbone. The PLC sequence generated in this way is then used 
to colour a homopolymer from the existing pool of good-solvent conformations, which is different from 
the one used to initiate the collapse process. This procedure is repeated 500 times, always ensuring 
that the homopolymer conformation that is decorated with the PLC sequence obtained at the end 
of the collapse and colouring process, is different from the starting homopolymer conformation. In 
this way, we ensure that we have a set of 500 PLC chains that differ from each other in sequence and 
in conformation. During the quench experiment, each of these PLC chains is then collapsed again, 
leading to 500 independent quench runs per parameter set. 

The second family of copolymers is the regular or alternating multi-block copolymer (MBC) with 
block length L (number of monomers in a contiguous H or P block). Hence the sequences were of 
the form (H^Pl^ where n = N/(2L) is the number of block pairs. No randomness is involved in 
this architecture; therefore again the sample size is approximately 500 independent quench runs per 
parameter set. 

The last family is known as the random-block copolymer (RBC) in which each block length L is 
an independent random variable sampled from a Poisson distribution 

/(L) = ^exp(-A), (9) 

where A = (L) is the average block length. The Poisson random variables were generated from a built- 
in function provided in MATLAB. For each of the <~ 500 stored starting conformations we generated 
one new sequence, such that again we have a total of <~ 500 independent quench runs per parameter 
set. We imposed a strict 1 : 1 constraint on the ratio of the number of monomers (H:P). Therefore, 
the H and P blocks were treated separately. Since the Poissonian random numbers usually do not add 
up to precisely N/2, we simply replaced the last random number with the remaining block. 
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2.3. Observables 

We monitor the time dependence of various observable quantities given below, which can be used to 
characterise the collapse kinetics. Important observables are the mean square radius of gyration 



and the internal energy 

U=J2 (Wr^c + VMrlu) [1 ~ ^c v }) , (11) 

where c M is unity for an H monomer and zero otherwise. 

The instantaneous shape of a polymer chain can be determined from its gyration tensor, G, whose 
Cartesian components are given by 

@ij = ^ t ( r A« — r cmi) ( r fj,j ~ r cmj) j (12) 
A" 

where the cm subscript refers to the chain's centre of mass. For each conformation, we determine the 
three eigenvalues A 2 , A 2 ,, A| of G, using the convention Af > A| > A§. The ratios of these numbers 
give an indication of the deviation of the molecule's shape from a sphere. 

The total collapse time r is defined as the time needed for the radius of gyration to reach 99% of 
its total change in size during the transition period in which the chain is transformed from the initial 
coil to the final state, i. e. 



R* s (T)-R* s>eq = —(R* s (0)-Rl eq ) (13) 



1 

fob' 

where R* cq = yJ(Rg 2 ) cq 1S the root mean square equilibrium dimensionless radius of gyration in the 
collapsed state (defined as the value obtained at the time at which the run was stopped). 
We also monitor the time dependence of the average cluster size, defined as 

(S n (t*)) = (14) 

using the algorithm of Ref. 48, where n(s) denotes the number of clusters of size s. Our previous 
investigation on homopolymers^ showed that this definition gives a more consistent picture of the 
dynamics than other possible definitions involving higher-order moments. Only H-type beads were 
used to define a cluster. Two H monomers are considered to be part of the same cluster if they (i) 
are not nearest neighbours along the backbone of the chain, and (ii) have an interparticle distance 
that does not exceed a certain value D, which we have set (in our units) to D* = 1.2 x 2 1 / 6 (a/l H ), 
following previous studies.^ 2 * 28 * 24 ! 

It is appropriate to note that all averages reported here are estimated by carrying out sample 
averages across the 500 independent runs, and the reported error is the estimated error in these 
averages. Except in the case of MBC chains, where the sequence is identical across all starting 
conformations, our data therefore involves double sampling over sequences and starting conformations. 
In order to avoid overcrowding in the plots, the statistical error bars are only included for a few selected 
points. In cases where the error bars are smaller than the symbol size, they are omitted. 
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5000 10000 15000 

t* 

00 

Figure 1. Variation of the mean square radius of gyration with time for N = 128 chains, in the 
presence of HI, for (a) the PLC chain and MBC chains with various values of block length L, and (b) 
all three types of copolymers with values of average block length (L) < 4. 
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3. Results and discussion 
3.1. Chain size 

Figure Q] shows the effects of the block length L for the protein-like, regular multi-block and random 
block copolymer chains on the evolution of the mean square radius of gyration for N — 128, in the 
presence of HI. Obviously, the average block length plays an important role in controlling the dynamics 
of the process as well as the final equilibrium size of the copolymer chains. For the PLC chains in our 
model, the average block length (L), for both N — 64 and N — 128 bead chains, is approximately 3, 
which is very close to the value reported for the analogous lattice modePM] ((L) = 3.173). 

We first discuss the effect of block length on the behaviour of the MBC chain. The simplest type 
of MBC is the diblock copolymer with L = N/2. It can be seen from Fig. Q] (a) that this diblock 
chain collapses to its equilibrium state faster than any of the other copolymer chains present in the 
plot. This is always true for such a copolymer because one of its ends (i. e. the end with the H block) 
behaves like a homopolymer chain in a poor solvent, which rapidly folds into its final equilibrium 
globular state, while the other end behaves like a homopolymer in a good solvent, remaining as a 
swollen coil. Due to the large size of this swollen part, the chain cannot form a compact equilibrium 
structure and consequently always has a relatively large final equilibrium size. 

By reducing the block size to L = 32, such that one has 4 blocks, one can see two time regimes of 
collapse. The first regime represents the early stage of collapse, where each H block of size L rapidly 
folds into a single cluster. Thus at the end of this stage, the chain consists of N j (2L) clusters or pearls 
separated by strings of P block chains. Since the folding of the H block is entirely homopolymer-like 
and is unaffected by the P type monomers, the time taken for each of these H block chains to completely 
fold into a small cluster is quite small, and it becomes systematically smaller for decreasing block size, 
such that it becomes hardly observable for shorter blocks. 

The second regime represents the growth or coalescence of these intrachain clusters into a single 
large cluster. These clusters are separated from each other by a string of P block, which has the same 
size as that of the original H block. If the separation between the two intrachain clusters is larger 
than the range of the attractive interaction (or, equivalently, if the P block length is large enough), 
then these intrachain clusters go through a diffusion-like process. It is only when they interact with 
each other, i. e. when their separation distance comes within the range of attractive interaction, that 
they coalesce. Thus the growth of the average cluster size is quite gradual in this case. However, if 
the P block chain length is sufficiently small, such that two intrachain clusters are within the range 
of their attractive interaction, then the two clusters quickly coalesce to form a larger cluster and 
diffusion plays little or no role. Moreover, the average cluster size will also grow much more rapidly, 
or, equivalently, there is a rapid reduction in chain size. Since the diffusion time increases with the 
size of the intervening P blockpSl it is expected that the coalescence time is shorter for smaller P block 
sizes and consequently the collapse is much more rapid compared to large P block sizes. Furthermore, 
a large P block size leads to a less compact final collapsed size because blocks of P beads are always 
repelled from the core of H beads, and they form dangling legs extruding away from the core. The 
results for L = 32 to L = 8 for MBC chains in Fig. Q] (a) also confirm these expectations. 

Note that for all MBC chains with L > 8, the final equilibrium size is still larger than that of the 
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Table 1. Values of the exponents for the early stages of collapse and the growth of the number 
average cluster size (i. e. a and z in (i?g(0)) — (i?g(i)) ~ t a and (S n ) ~ t z ) for all types of copolymers 
with N — 128 at various values of average block lengths (L), in the presence of HI. Values of the total 
collapse time r and the equilibrium mean square radius of gyration (-Rg) are listed as well. Note that 
PLC, MBC and RBC denote protein-like, multi-block and random-block copolymers, respectively. 

PLC chain. However, when L — 4, the MBC chain reaches its equilibrium native state faster, and 
the final equilibrium size is smaller than that of a PLC chain. Reducing the block size of the MBC 
chain to L = 3 further reduces the final size as well as the total collapse time. In fact, the data shows 
that the MBC chain with L — 3 produces the lowest final equilibrium size compared to all other types 
of copolymer chains for the entire range of L and (L) values that have been investigated. The total 
collapse time and the mean square radius of gyration for all three types of copolymers are listed in 
Tabled! 

Further reduction in the block size to L = 2 leads to a slightly faster collapse (as evidenced by 
the value of r in Table [T]), but it also increases the final chain equilibrium size. Note that for every 
L monomers of H type that fuse with another H cluster, L monomers of P type are brought into 
close proximity, due to the connectivity along the chain. As a result, a repulsive force builds up as 
the coalescence process takes place. For the L = 2 chain, the energy gained from coalescence cannot 
overcome this repulsion, and therefore the chain does not form a single cluster of H monomers, i. e. 
it cannot fold into a spherical compact globule. The L — 1 case conforms to this behaviour, with the 
chain collapsing slower, and the final equilibrium size being even larger than for L = 2. 

The large repulsive force built up due to the overcrowded presence of P type monomers for MBC 
chains with L < 2 can also be observed via the high value of the internal energy U in Fig. [2j Note 
that the internal energy for MBC chains with L = 32 and 16 is not reported here because the chains 
with these block sizes had not yet reached their equilibrium state, which is also subject to large 
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Figure 2. Coordinate pairs (U, {R g 2 } ) of the internal energy and (R g 2 ) , respectively, in the final 
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collapsed state for all copolymer chains with N = 128, in the presence of HI. 



fluctuations. This figure also clearly indicates that a chain which has the lowest equilibrium size does 
not necessarily have the lowest internal energy. Moreover, a chain with the lowest internal energy may 
not fold into the most compact final structure. This finding is consistent with the results of Ref. l28l 
where BD simulations without HI were performed. 

From Fig. [T] (b), we can see that the PLC chain and the RBC chain with (L) — 4 are almost 
identical in terms of their final sizes and their time behaviour of R g . However, reducing the average 
block size of the RBC chain to (L) = 3 and (L) = 2 systematically decreases both the collapse time 
and the equilibrium size even further. This finding is at variance with the results of Ref. 1371 where 
the PLC chain was found to collapse more rapidly and into a more compact structure than the RBC 
chain with the same (L). It is not quite clear what the origin of this discrepancy is; most likely it 
is the fact that our model is a bead-spring model in the continuum, while Ref. [37] studied a lattice 
model, which results in different local packing geometries, which are certainly very important for this 
problem. Another possible source could be different interaction parameters. The fact that we study a 
system with HI and Ref. |37jone without can however not be used as an explanation; as Fig. [3] shows, 
we find the same behaviour also in the absence of HI (see below). For the protein folding problem 
this result indicates that the relation between the sequence and the thus-encoded collapse kinetics 
and native state is probably more intricate than the simple and appealing picture that the labeling 
procedure of Ref. l37l suggests. 

We fitted a power law 

(Rl(t)) = (Rl(0))-At a (15) 
to the initial decay of the mean square gyration radius, and the resulting exponents a are listed 
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Figure 3. Variation of the mean square radius of gyration with time for all three types of copolymer 
chains with N = 128, for various values of the average block length (L), in the absence of HI. 

in Table [T] For a homopolymer chain at the same quench depth, Ref. [23] had obtained a value of 
a = 1.05 ± 0.01 in the presence of HI. Interestingly, it can be seen from Tabled] that there exist some 
copolymer chains which have a larger exponent for this early stage compared to that of a homopolymer 
chain. Similar results have also been observed in Ref. [55] for copolymer chains in the absence of HI. 

3.2. Influence of HI 

In order to study the influence of HI, we have also carried out simulations without HI, for the case N = 
128, and the results are shown in Fig. [3] We find that in general the collapse dynamics is substantially 
slowed down, compared to the HI case, which corresponds to the observation for homopolymers that 
the presence of HI (Zimm model) leads to systematically faster dynamics than in the case when HI 
is absent (Rouse model). In most cases, the chains had not even relaxed fully into equilibrium at 
the time at which the HI runs (and also the non-HI runs) were stopped. For the one case in which 
equilibrium was reached, we find the same value for (-R|) as in the case with HI, as it must be. The 
qualitative results concerning the influence of the chain architecture on the collapse kinetics remain 
however completely unchanged: Again, the PLC chain corresponds to the (L) = 4 RBC, but the 
RBCs with shorter block length collapse even faster. 

The observed acceleration of the collapse dynamics as a result of HI is consistent with our previous 
findings on homopolymer collapse pi where it was shown that for N — 128 the inclusion of HI reduces 
the total collapse time by at least a factor of two. A complete comparison between the total collapse 
times for cases with and without HI is not reported in this work, due to the slow convergence to 
equilibrium in the absence of HI. The only case where equilibrium was reached (MBC, L — 2) seems 
to indicate that the speedup is even larger in the present case (r = 5233 without HI, r = 1967 




12 



^^^^ 




> PLC, e*=2.0 
PLC, e*=2.5 

o PLC, e*=5.0 

+ PLC,e*=10.0 

* MBC, e*=2.0 
a MBC, e*=2.5 

* MBC, e*=5.0 

* MBC,e*=10.0 



2000 



4000 

t* 



6000 8000 10000 



(a) 




C> PLC, e*=2.0 
PLC, e*=2.5 

o PLC, e*=5.0 

+ PLC,e*=10.0 
RBC, e*=2.0 

□ RBC, e*=2.5 
RBC, e*=5.0 

x RBC,e*=10.0 



2000 4000 6000 8000 10000 

t* 



(b) 



Figure 4. Variation of the mean square radius of gyration with time for chains with N = 64, in the 
presence of HI, for (a) PLC and MBC chains with block length of L — 3, and (b) PLC and RBC 
chains with an average block length of (L) = 3. Here, e* denotes the value of the quench depth or 
^hj/k^T. Insets: Variation of the mean square radius of gyration with time for chains with N = 64 
and L = 3 or (L) = 3 at various quench depths e*, for (a) MBC chain, and (b) PLC chain. 



13 

with HI). Since hydrodynamics is a large-scale phenomenon, one expects the effect of HI to be the 
more dramatic the more extended the objects are. Due to the uncollapsed P strings, the chains stay 
systematically longer in an extended state than in the homopolymer case, and this is probably the 
reason why HI has an even stronger effect. 

3.3. Influence of quench depth 

We have simulated N = 64 chains at various values of the quench depth e* , in order to study the 
influence of this parameter, while keeping the block length (or average block length) at L — 3 (or 
(L) = 3). Figure S] shows the results for the gyration radius. For all the quench depths investigated, 
it can be seen that the MBC chains collapse the fastest and have the most compact equilibrium size, 
followed by the RBC chains. Out of all three types of copolymer chains, the PLC chain collapses 
the slowest and the final equilibrium size is the least compact. Generally, one would expect that 
the radius of gyration for a chain at deep quench should be smaller than that of a chain with lower 
quench values because the stronger attractive interaction would cause the chain to squeeze into a 
tighter globule. The above figure shows that a chain smoothly folds into its final compact state for 
low quench depth. However, at very large quench depth (i. e. t^j/k^T = 10), a chain gets trapped 
in a metastable state and stays there for a long time rather than approach its final minimum energy 
state. The inset of Fig. 3] (a) shows this trapping behaviour at very large quench depths much more 
clearly for an MBC chain. This inset reveals that, while MBC chains with e^/k-oT < 5 have fully 
reached their equilibrium compact state, chains with t^j/k^T = 10 are still gradually approaching the 
equilibrium state. A similar result is also observed for other types of copolymer chains (see Fig.@](b)). 
A similar behaviour has been observed by various other authorsl±mil2iii2| i n R e f. [M]it was found that 
a homopolymer chain will get trapped at t^j/k^T — 5, while our results indicate that a copolymer 
chain still smoothly folds at this quench depth. This result seems to indicate that the presence of P 
type monomers in the chain prevents it from being trapped in a local well and smooths out the energy 
landscape for the folding process of copolymers. This in turn pushes the value of quench depth where 
trapping occurs to a much higher value. 

3.4. Snapshots 

Snapshots of the typical kinetic pathways for a collapsing chain of length N — 128 for MBC chains are 
shown in Fig. [5l and for PLC and RBC chains in Fig. [5] For some selected cases, we show the time 
evolution in the movie files Sl.avi and S2.avi of the Supporting Information. It can be seen that there 
exist at least two distinct stages for the kinetics of collapse for these copolymers. It is to be noted that 
almost all of the copolymer chains with sequences that fold into a spherical compact structure have 
a three-stage mechanism. The early stage involves the rapid formation of localised globules along 
the chain. Following this stage, these globules then coalesce together by fusing with other nearby 
blobs to form a dumbbell or two pearls separated by linear chain. Finally, the pearls combine to form 
a sausage which slowly rearranges itself into a compact state. However, for copolymer chains with 
sequences consisting of short P blocks that do not fold into a spherical compact structure (for instance, 
MBC chains with L < 2), only the first two distinct stages are observed and the final compaction 
stage apparently does not occur. Furthermore, copolymer chains with sequences consisting of large 
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Figure 5. Snapshots of different types of collapsing regular multi-block copolymers (MBC) chain 
with N — 128, in the presence of HI. 
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Figure 6. Snapshots of different types of collapsing protein-like copolymers (PLC) and random-block 
copolymers (RBC) chains with N = 128, in the presence of HI. 




t* 



Figure 7. Variation of the ratio of the largest eigenvalue to the smallest eigenvalue of the gyration 
tensor with time for all three types of copolymer chains with N = 128, for various values of the average 
block length (L), in the presence of HI. Inset: The corresponding variation of the smallest eigenvalue 
with time. Here, the symbol A, is actually an abbreviation for (A^) 1 2 . 
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P blocks, such as MBC chains with L = 32 and 16, seem to acquire another distinct stage after the 
early rapid collapse stage, known as the "diffusion" or "plateau" regime. During this diffusion stage, 
the number and size of the intrachain clusters remains unchanged due to the large separation between 
the H blocks. Similar qualitative features of the collapse pathways have also been observed by other 
authors™ 

3.5. Asphericity 

Apart from the snapshots, the asphericity of the chain can also be observed from the eigenvalues of 
the gyration tensor. Figure [JJ shows the time evolution of the ratio of the largest eigenvalue to the 
smallest eigenvalue for all three types of copolymer chains with N — 128, for various values of the 
average block length (L), in the presence of HI. Interestingly, the data in Fig. [JJ clearly shows that this 
ratio of eigenvalues behaves non-monotonically with time. A similar result is also observed for the 
ratio of the intermediate eigenvalue to the smallest eigenvalue A2/A3. The non-monotonic behaviour 
is mainly due to the non-monotonic nature of the smallest eigenvalue A3 , as can be seen in the inset 
of Fig. [JJ This behaviour can be explained by the fact that during the collapse of the chain, there is 
a sudden change in the shape of the chain from an initial "ellipsoid" (note that a typical self-avoiding 
walk is highly anisotropic) to a pearl-necklace and finally (typically) to a sphere. During this shape 
transformation, the eigenvalue for the smallest principal axis is first reduced to the lowest value, at 
which the thinnest pearl- necklace chain is formed. This principal axis then grows in width when the 
pearl-necklace chain starts to transform by absorbing nearby pearls or clusters, and when the bridging 
strings of monomers form a final spherical structure. Note that a similar result has also been observed 
for a homopolymer chainPH 

3.6. Cluster analysis 

The growth of the number average cluster size with time for our chosen value of overlapping distance 
D* = 1.2 x 2 1 ' 6 {a/lff) is shown in Fig. [5J The most striking aspect about these data is their quanti- 
tative demonstration of the multi-stage process of collapse kinetics, already indicated qualitatively by 
inspection of snapshots, in a much clearer fashion than by the data on the gyration radius. Typically, 
there are three growth regimes: (i) Collapse of the single H blocks on a fairly short time scale, giving 
rise to a moderate increase in the mean cluster size (S n ); (ii) Coalescence of these clusters at a later 
time scale, resulting in a steep increase in cluster size; and finally (iii) internal rearrangement of the 
resulting "sausage" that does not change (S n } anymore. However, for MBC chains with large block 
size (i. e. L = 32 and 16, see Fig. [5] (c)) yet another (fourth) regime can be observed, which follows 
directly the initial single-block collapse, and which is characterised by a constant value of (S n ) . This 
clearly corresponds to the time regime that is needed for the single globules to find each other in a 
diffusion-like process. 

This diffusion regime for MBC chains with large P blocks is also observed in the time evolution 
of the number of clusters N c , as shown in Fig. |H1 While the L ~ 64 chain (diblock copolymer) shows 
a simple decay of iV c to N c = 1, i. e. coalescence of smaller clusters to one large globule, the chains 
with more blocks (L = 32, 16, 8) exhibit a clear plateau at N c = N/ (2L), i. e. the number of H blocks. 
For smaller block sizes, such a plateau is no longer clearly identifiable, since the intervening P blocks 
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Figure 8. Variation of the number average cluster size with time for chains with N = 128 in the 
presence of HI, for (a) MBC chains, (b) the PLC chain and RBC chains with (L) < 4, and (c) MBC 
chains with L = 32 and 16. 
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Figure 9. Variation of the number of clusters with time for the PLC and MBC chains, in the presence 
of HI. 




Figure 10. Coordinate pairs (Lp, z) of the block size of P type monomers and the coarsening stage 
exponent, respectively, for MBC chains and a homopolymer chain with N = 128, in the presence of 
HI. 
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are too short to keep the H clusters away from each other. Interestingly, for L = 2, N c exhibits a 
clear maximum, which corresponds to breakup and coalescence of clusters. Furthermore, the PLC 
chain, whose data is shown for comparison, shows just a smooth decay to iV c = 1, without any visible 
structure. 

During the cluster coalescence stage, previous studies of homopolymer chains indicate that one 
might expect a power law 

(S n (t))~t z . (16) 

We therefore fitted such a behaviour to our data; the results for the exponent z are listed in Table [TJ 
Interestingly, the z values that we find for our copolymers are always smaller than the value z = 
1.08 ± 0.01 obtained for a homopolymer P2 Given the fact that the obtained values differ from each 
other substantially, it is far from clear that the obtained numbers have much to do with a well-defined 
asymptotic power law. Furthermore, the fact that the coalescence process involves two competing 
length scales (the average globule size of the consolidated clusters, and the size of the P "loops") makes 
the application of self-similarity arguments, which are usually employed in order to justify power-law 
behaviour, rather doubtful. On the other hand, the fact that the z value depends systematically 
on the block length (as shown in Fig. [TU] for the homopolymer and MBC chains up to block length 
L = 32) makes the interpretation of z in terms of an effective exponent, which in reality describes a 
more complicated dynamics, fairly likely. 

Noticeably, data for MBC chains with L — 1 and L = 2 in Fig. [10] do not lie on the smooth curve 
that can be drawn by eye through all the remaining chains starting from the homopolymer chain to 
the MBC chain with L = 32. As discussed earlier, for these two MBC chains, repulsive forces build up 
during the coalescence process due to the close proximity of the P blocks with each other. This leads to 
an equilibrium state with the highest energy (Fig. [5]), and the highest asphericity, as evident from the 
snapshot (Fig. [5]), and the ratio of the largest to smallest eigenvalue (Fig. [T])- The non-occurrence of 
the final cluster consolidation stage into a compact globule as a consequence of the built up repulsive 
energy is probably the reason for the data for these two chains not following the general trend seen 
in the case of the other MBC chains. A similar pattern is seen below in the correlation between the 
collapse time and block size, and equilibrium mean square gyration radius. 

3.7. Collapse time 

The characteristic collapse times (r) for all three types of copolymers with N = 128 are shown in 
Table [TJ It is observed that for all types of copolymers used in this work, the collapse time is much 
larger than the collapse time for a homopolymer chain, for which r = 591 ± 7. As has been pointed 
out earlier, MBC chains with L = 32 and 16 had not yet reached their equilibrium state during the 
course of our simulation, and hence the total collapse time for these two chains are not available. 
The collapse time also increases systematically with the block size, as demonstrated in Fig. [TT] for the 
homopolymer and MBC chains up to L = 8. These data hence clearly indicate that the block size of 
P monomers not only controls the duration of the diffusion process, but also the collapse rate of the 
cluster coalescence stage for copolymers with small P block size where diffusion plays little or no role. 
The combined diffusion and cluster consolidation process is thus identified as the decisive rate-limiting 
factor for the overall collapse. The absence of these processes for a homopolymer leads to a much 
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Figure 11. Coordinate pairs (Lp, r) of the block size of P type monomers and total collapse time for 
MBC chains and a homopolymer chain with N = 128, in the presence of HI. 
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Figure 12. Total collapse time versus equilibrium mean square gyration radius for N — 128 copolymer 
chains of different types, in the presence of HI. 



21 




Figure 13. The tails of the normalised distribution of P blocks for a protein-like copolymer (PLC) 
and random-block copolymers (RBC) with three different average block lengths (L), for N = 128. 
Inset: The complete normalised distribution for PLC and RBC chains. 

faster collapse, and also leads in the case of MBC chains with L = 1 and L = 2, to a departure from 
the observed trend for the other chains. 

In order to understand the relationship between the kinetic accessibility of the final state and the 
final equilibrium size, we plot the collapse time versus the final size in Fig. [T^J Except for MBC chains 
with L = 1 and L = 2 (for the reasons discussed earlier), the figure clearly demonstrates that the total 
collapse time is directly related to the final equilibrium size. This result reveals a very interesting 
feature that is somewhat unexpected, i. e., a chain with a small equilibrium size tends to fold much 
more rapidly than a chain with a larger equilibrium size. Intuitively, one might expect that it would 
take longer for a chain to fold into a more compact equilibrium structure rather than into a loosely 
packed structure, but the present results suggest the opposite. In Ref. [51] it was pointed out that a 
pronounced energy minimum is a necessary condition to guarantee that the native state is stable, and 
that such a minimum is sufficient, for a compact globule with random structures, to rapidly find the 
native state. Thus the deep energy minimum of the compact structure seems to provide a guide (or a 
strong thermodynamic driving force) for the chain to quickly fold into its final equilibrium state. In 
the present work, we have made no attempts to determine the free energy of the chains. However, we 
anticipate that knowledge of the relationship between the free energy of the native state and the size 
of the native conformation may help in understanding this behaviour better. 

3.8. Block size distribution 

As has been discussed above, the length of P blocks for MBC chains has a very decisive influence on the 
kinetics (large blocks lead to a substantial slowing down of the process) as well as on the equilibrium 
size (large blocks increase the size). Since this influence is very strong, one might expect that in 
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disordered copolymers the presence of such large blocks should also have a substantial influence, 
even if their probability of occurrence is fairly small. We have therefore analysed the probability 
distribution of P block lengths for the PLC chains that we have generated, and compare it in Fig. [T3J 
to the Poisson distributions of RBC chains with similar (L). The inset shows the total distributions, 
but for our purposes the tails of the distributions at large L values are of particular interest; hence 
they are shown in the main part. We find the very interesting result that for the PLC chains the tail 
is substantially enhanced compared to the Poissonians with 2 < (L) < 4. It is therefore quite natural 
to assume that it is this enhanced tail that leads to the slowed-down dynamics of the PLC, compared 
to the RBC with identical (L). If this correlation is really valid, then one should expect that the 
lattice PLC chains of Ref. [37] should show a deflated tail, in comparison to the corresponding RBCs, 
as a result of different packing behaviour. It should be quite interesting to test this. 



4. Summary and conclusions 

We have shown that Brownian Dynamics simulations incorporating implicit hydrodynamic interactions 
can be used to study the dynamics of copolymer collapse in a poor solvent. Our observations regarding 
the speedup of collapse caused by hydrodynamic interactions, and the existence of at least two stages 
of collapse, are similar to those that have been reported previously in the literature. Due to the 
inherent computational advantages of an implicit-solvent model, compared to, say, straightforward 
Molecular Dynamics with explicit solvent, we were able to study the phenomena with reasonable 
statistical accuracy. 

The kinetics of collapse can be described as a rapid initial formation of clusters followed by cluster 
coalescence and sometimes a rearrangement of the final cluster to form a compact state. It is also found 
that the presence of P monomers pushes the value of quench depth at which trapping phenomena occur 
to a higher value compared to the value seen for homopolymer chains. A striking feature observed 
here is that the total collapse time is completely governed by the cluster consolidation stage and the 
rate of collapse of this stage depends on the block size of P monomers along the chain. It is found that 
in our model random block copolymers collapse more rapidly than "protein-like" copolymers with the 
same average block length, if the latter are generated by the labelling procedure suggested in Ref. [37J 
and this behaviour can apparently be traced back to an enhanced occurrence of long P blocks in the 
PLCs. 

We finally emphasise again that our model is intended as a "toy model" designed to study the 
influence of certain aspects (monomer sequence, hydrodynamic interactions) on the collapse kinetics, 
which cannot be expected to be directly related to real proteins. To illustrate this point, let us assume 
that each model monomer corresponds to a single amino acid residue, with a typical size of, say, Iran, 
which would then be the value of Ih in real units. Since we are only interested in order-of-magnitudc 
estimates here, this could also be taken as the value for a or for Qo. A chain of N — 64 or 128 would 
then correspond to a fairly short polypeptide whose native state would typically need to be described 
in terms of a much more elaborate geometry (e. g. alpha helices, beta sheets) than just a spherical or 
sausage-like globule. The elementary time scale Xh is then found via 

A * = 4^ 6 ^|r' < 17 > 
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where a, the Stokes radius of the model monomer, is given by 

a = h* 'y/nljj, (18) 

such that we find 

Xh = I^St- (19) 

For our value of h* , and assuming the viscosity of water and room temperature, we thus find A# ~ 
0.3ns as our elementary time unit. A collapse time of 6000 in our BD units would thus correspond to 
roughly 2/is. This value should however not be taken too seriously, since our model does not include 
bond-bending and torsional potentials and other interactions that would be needed to describe the 
structure and energetics on such a small length scale. Furthermore, the value of Ih, which we have 
essentially just guessed, enters with the third power. Nevertheless, it may be noted that detailed 
simulations and experiments on real proteins^ have also produced collapse times in the microsecond 
range, whose values are however somewhat larger, for chains that are even shorter. In general, we 
believe that the simplicity of our model results in a lack of local kinetic barriers, such that our estimate 
for the collapse time tends to be smaller than the folding time of the "corresponding" real protein. In 
this context, it should also be noted that the interpretation of one model monomer in terms of just 
one residue is certainly not imperative. One could as well consider a model where one model monomer 
corresponds to several residues (under the assumption that proteins with such long blocks exist). In 
this case, one would have a larger value for Ih, and therefore a very much larger folding time in real 
units. 
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